# load and attach packages----
library(ggplot2)
library(margins)
library(plm)
library(haven)
library(jtools)
library(Zelig) # Can be used for rare events logit, but cannot handle robust clustered errors :(
library(bife) # For fixed effects logit
library(pcse)
library(visreg) # Plots Interactions and predicted values against residuals
library(effects) # Plots effects including interactions
library(tidyverse) # For data manipulation
library(dplyr) # For data manipulation
library(Hmisc) # for summary statistics
library(foreign) # allows the exporting of R back to Stata

# These packages are for making MSWord tables
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(devtools)
library(PanelMatch)
library("car")
library(ggplot2)
library(tidyr)
library(plyr)

# Load Data----Enter the full path for local location of data
SaCrat_data01.2 <- read_dta("SanctionsAndStockMarkets.dta")

# convert SD and countyrmo to integers. findAllTreated needs them to be integers later
SaCrat_data01.2$ccode <- as.integer(SaCrat_data01.2$ccode)
SaCrat_data01.2$countyrmo <- as.integer(SaCrat_data01.2$countyrmo)
SaCrat_data01.2<-subset(SaCrat_data01.2, (!is.na(SaCrat_data01.2$ImpD1)))
SaCrat_data01.2 <- data.frame(SaCrat_data01.2)

#PanelMatch 
# With CBPSW refinement
PM.resultsI.CBPSW <- PanelMatch(lag = 12, time.id = "countyrmo", unit.id = "ccode", 
                          treatment = "ImpDT5", refinement.method = "CBPS.weight", 
                          data = SaCrat_data01.2, match.missing = T, 
                          covs.formula = ~ I(lag(StMk, 1:3)) + IntCon +CW +prelec +ideoLef +openk +EXRat +CPI +lnTR +lnGDP
                          + NotAnyIMPTty + Dev_Cou, 
                          size.match = 10, qoi = "att" ,outcome.var = "StMk", verbose=TRUE, 
                          lead = 0:3, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

# With no refinement
PM.resultsI.none <- PanelMatch(lag = 12, time.id = "countyrmo", unit.id = "ccode", 
                         treatment = "ImpDT5", refinement.method = "none", 
                         data = SaCrat_data01.2, match.missing = T, 
                         covs.formula = ~ I(lag(StMk, 1:3)) + IntCon +CW +prelec +ideoLef +openk +EXRat +CPI +lnTR +lnGDP
                         + NotAnyIMPTty + Dev_Cou, 
                         size.match = 10, qoi = "att" ,outcome.var = "StMk", verbose=TRUE, 
                          lead = 0:3, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

# generate panel estimates With CBPS.weight Refinement
PE.results95I.CBPSW <- PanelEstimate(inference = "bootstrap", sets = PM.resultsI.CBPSW, 
                                     data = SaCrat_data01.2, CI = .95, ITER = 1000, set.seed(604))
summary(PE.results95I.CBPSW)

PE.results90I.CBPSW <- PanelEstimate(inference = "bootstrap", sets = PM.resultsI.CBPSW, 
                                     data = SaCrat_data01.2, CI = .90, ITER = 1000, set.seed(605))
summary(PE.results90I.CBPSW)

#################################
# Plot ATT with 90% and 95% standard errors -- CBPSW
#################################
results95I.CBPSW <- as.data.frame(summary(PE.results95I.CBPSW)$summary)
results90I.CBPSW <- as.data.frame(summary(PE.results90I.CBPSW)$summary)

theme_set(theme_bw(base_size = 15))

# PLOT THREE MONTH WINDOW
I <- ggplot(data = results95I.CBPSW, aes(y = estimate, x = 1:4)) +
  geom_errorbar(aes(ymin = `2.5%`, 
                    ymax = `97.5%`), width=0, size = 1) +
  geom_errorbar(data = results90I.CBPSW, aes(ymin = `5%`,
                                             ymax = `95%`), 
                width=0, size = 2) +
  geom_point(size = 6.) +
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.5) +
  scale_y_continuous() +
  scale_x_discrete(name="Months Since Import Sanction", limits=c("0"="t+0","1"="t+1","2"="t+2","3"="t+3"))+
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank()) +
  ylab("Estimated Effect of Treatment on the Treated")

I

#################################
# CHECK AND PLOT CHANGE IN COVARIATE BALANCE
#################################
# Calculate covariate balance
refined_bal <- as.data.frame(get_covariate_balance(PM.resultsI.none$att, SaCrat_data01.2, 
                                                   covariates = c("IntCon", "CW", "prelec", "ideoLef", "openk",
                                                                   "EXRat", "CPI", "lnTR", "lnGDP", "NotAnyIMPTty",
                                                                  "Dev_Cou" ), 
                                                   plot = F))
# Reformat for plotting
refined_bal <- gather(refined_bal, key = "covariate", value = "ref_sd")
refined_bal$t <- 12:0
refined_bal$abs_ref_sd <- abs(refined_bal$ref_sd)

# Same procedure for "unrefined" sets
unrefined_bal <- as.data.frame(get_covariate_balance(PM.resultsI.CBPSW$att, SaCrat_data01.2, 
                                                     covariates = c("IntCon", "CW", "prelec", "ideoLef", "openk",
                                                                     "EXRat", "CPI", "lnTR", "lnGDP", "NotAnyIMPTty",
                                                                    "Dev_Cou"), 
                                                     plot = F))
# Reformat for plotting
unrefined_bal <- gather(unrefined_bal, key = "covariate", value = "unref_sd")
unrefined_bal$t <- 12:0
unrefined_bal$abs_unref_sd <- abs(unrefined_bal$unref_sd)

# Merge into a single dataframe
balance_test <- merge(refined_bal, unrefined_bal, by = c("covariate", "t"))
balance_test$covariate <- as.factor(balance_test$covariate)

# Create a reference line
refLine <- data.frame(x = c(0, 1), y = c(0, 1))

# Plot the covariate balance of "refined" set against "unrefined" set
ggplot(data = balance_test, 
       aes(x = abs_ref_sd, y = abs_unref_sd, color = covariate)) +
  geom_point(size = 4) +
  scale_y_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_line(data = refLine, aes(x = x, y = y, color = NULL), linetype = 2) +
  ylab("Post-Refinement") +
  xlab("Pre-Refinement")

rm(refLine, refined_bal, unrefined_bal, balance_test)
